#include "covMatrix.h"

/**
* @brief 矩法求归一化协方差矩阵
* @param infile  				译csv格式文件保存的原始数据，列数为nDim 
* @param n   					输入数据列数 
* @param covMatrix       		长度为nDim*nDim的数组，存放返回的协方差矩阵
*/
void c_covMatrix(FILE *infile, int n, double *covMatrix)
{ 
  double *row=malloc(n*sizeof(double));
  double *s1=malloc(n*sizeof(double));
  double *s2=covMatrix;
  memset(s1,0,n*sizeof(double));
  memset(s2,0,n*n*sizeof(double));
  int cnt=0;
  
  int bufLen=n*n*5*sizeof(char);
  char *buf=malloc(bufLen);
    
  while (1) 
  { 
    buf=fgets(buf,bufLen,infile);  
    if (feof(infile)) break;
    
    char *h=buf;
    char *p=buf;
    double *r=row;
    // 拆分行转成double数组
    while (*p != '\0')
    { 
      if (*p != ',') 
        p++; 
      else
      { 
        *p='\0';
        *r=atof(h);
        p++;
        h=p;    
        r++;
      }
    }
    *r=atof(h);
    
    // 累加count，sum(x),sum(xi*xj) 上三角矩阵
    cnt++;
    for (int i=0; i<n; i++)
    { 
      s1[i] += row[i];
      for (int j=i; j<n; j++)
        s2[i*n+j] += row[i]*row[j];
    }
  }
  
  // E[x]和E[xi*xj]
  for (int i=0; i<n; i++)
  { 
    s1[i] /= cnt;
    for (int j=i; j<n; j++)
      s2[i*n+j] /= cnt;
  }
  
  // 协方差E[xi*xj]-E[xi]*E[xj]
  for (int i=0; i<n; i++)
    for (int j=i; j<n; j++)
      s2[i*n+j] -= s1[i]*s1[j];
  
  // 归一化协方差S[xi*xj]/srqt(S[xi]*S[xj])
  for (int i=0; i<n; i++)
    s1[i] = sqrt(s2[i*n+i]);

  for (int i=0; i<n; i++)
    for (int j=i; j<n; j++)
    { 
      s2[i*n+j] /= s1[i]*s1[j];
      s2[j*n+i] = s2[i*n+j]; //协方差是对称矩阵
    }
    
  free(row);
  free(s1);
  free(buf); 

  return;
}
